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ABSTRACT 

An  analysis  is  made  of  the  analogy  between  schemes  for  integrating 
the  equations  of  motion  of  a  compressible  fluid,  in  Lagrangian  coordi¬ 
nates,  and  molecular  models  of  matter.  Computation  schemes  involving 
one  and  two  independent  variables  are  considered,  for  flows  with  and 
without  shocks. 
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V 


I.  FLOWS  IN  ONE  SPACE  DIMENSION 

1 •  Introduction.  In  his  pioneer  paper  [9]  on  the  numerical  solu¬ 
tion  of  time -dependent  compressible  flow  problems,  von  Neumann  suggested 
using  Lagrangian  independent  variables.  That  is,  he  suggested  computing 
the  vector  position  x  of  a  moving  fluid  particle  with  material  coordinate 
a,  as  a  function  x(a,t)  of  time  t.  For  homentropic  flow  in  one  space  di¬ 
mension  (see  §  2),  this  amounts  to  integrating  the  hyperbolic  equation 

(1)  Xtt  =  -  pa  =  F'(a)  xaa'  where  a  =  xa; 
and 

(2)  P  =  Pc  -  F(o) 

is  a  homentropic  equation  of  state.  Here  subscripts  denote  partial  deriv- 

2  2  1 

atives  (thus  x^  =  d  x/dt  ),  a  =  p  is  the  specific  volume,  and  F’(o)  >  0. 

* 

Von  Neumann  observed  that  the  natural  spatial  discretization  of  (l) 
yielded  a  system  of  ordinary  differential  equations  of  the  form 

(3)  m  d^/dt2  =  f(Xj.  -  xi_1)  -  f(xi+1  -  x±), 


Though  [9]  is  not  easily  available,  it  has  been  summarized  in  [6], 


V 
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where  we  set  f(5)  =  -F(&/Aa).  These  are  the  equations  of  motion 

of  a  system  of  identical  particles  of  mass  m  =  A  a,  joined  by  springs  — 

i.e.,  of  a  mass -spring  system. 

Von  Neumann  also  suggested  that  the  trajectories  x^t)  of  these 
particles  were  analogous  to  those  of  gas  molecules  in  the  kinetic  theory 
of  gases.  Essentially  because  of  this  analogy,  he  "ejected  that  the 
approximation  (3)  would  prove  itself  better  than  its  original  (l),  and 
give  adequate  approximate  descriptions  of  shocks". 

Since  1 9k-b,  von  Neumann's  computational  ideas  have  been  extensively 
explored.  Von  Neumann  himself  [11]  decided  to  introduce  "artificial  vis¬ 
cosity"  (dashpots)  into  the  mass-spring  system  (3)  to  simulate  shocks, 
as  being  more  efficient  computationally  than  statistical  averaging  of 
the  type  used  in  the  kinetic  of  gases.  But  no  systematic  exploration  has 
been  made  of  the  analogy  which  seemed  so  interesting  to  von  Neumann. 

The  present  paper  represents  such  an  exploration;  it  analyzes  the 
relation  between  Lagrangian  hydrodynamical  computation  schemes  and 
molecular  models  of  matter. 

First,  it  notes  that  the  force-laws  of  the  system  (3)  correspond 
to  molecular  models  of  crystals  and  liquids  ("dense"  gases)  more  closely 
than  to  the  models  of  "rare"  gases  used  in  the  kinetic  theory  of  gases. 

It  then  shows  how  various  physical  concepts,  familiar  in  theories  of  the 

*It  was  begun  jointly  with  S.  Ulam  in  1953,  in  informal  collaboration 
with  von  Neumann.  The  collaboration  of  Mr.  Lynch  started  in  1958* 

The  work  has  been  done  largely  at  the  Los  Alamos  Scientific  Laboratory , 
and  trnder  Contract  AT(30-1  )-1987* 
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liquid  and  solid  state,  are  applicable  to  hydrodynamical  computation 
schemes.  It  continues  with  applications  of  ideas  from  statistical  me¬ 
chanics  to  such  computation  schemes. 

Finally,  it  explores  the  mechanical  properties  which  are  possessed 

t>y  the  synthetic  materials  defined  by  such  computation  schemes.  Such 
* 

materials  may  well  simulate  the  behavior  of  real  solids  under  large  de¬ 
formations,  and  that  of  real  solids  and  liquids  under  extremely  high 
rates  of  strain,  as  well  as  or  better  than  the  classical  models  of  con- 

y  y 

tinuum  mechanics.  In  addition,  their  properties  can  be  effectively  in¬ 
vestigated  on  high-speed  computing  machines,  and  a  number  of  specific  in¬ 
vestigations  of  this  type  are  proposed. 

2.  Elastic  fluids.  The  simplest  model  of  continuum  mechanics  is 
defined  by  the  following  two  assumptions:  (i)  the  stress  at  any  point 
x  and  time  t  is  a  scalar  pressure  p(x,t),  and  (ii)  the  stress  is  a 
single-valued  function  of  the  density,  p  =  p(p).  These  assumptions  may 
be  said  to  define  a  homogeneous  elastic  fluid,  since  they  hold  in  any 
homogeneous  isotropic  elastic  medium  incapable  of  sustaining  a  static 
shear  stress. 

In  fluid  mechanics,  flows  with  these  properties  are  called  homen- 
tropic.  They  arise  in  the  isentropic  (e.g.,  adiabatic)  flow  of  a  homo¬ 
geneous  fluid.  We  omit  the  proof,  and  also  the  derivation  of 

In  the  limiting  case  At  -* 0  of  zero  time-step,  they  constitute  molecular 
models  of  matter  in  the  classical  sense. 

** 

To  simulate  quantum-mechanical  phenomena,  however,  one  must  use  tech¬ 
niques  very  different  from  those  discussed  below. 
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Eqs.  (l)-(2)  of  §  1,  as  a  convenient  mathematical  expression  of  these 
hypotheses. 

In  the  adiabatic  flow  of  a  perfect  gas,  we  have  PQ  =  0  and 
(4)  p  =  -  F(o)  =  kef7,  F'(a)  =  kyO~7~  , 

the  so-called  "polytropic"  equation  of  state.  (See  the  discussion  in  §  4.) 
The  approximation  (4)  seems  physically  adequate  for  atmospheric  air,  with 
y  =  1.4.  A  wide  variety  of  other  gases  can  also  be  fitted  quite  closely 
by  (4) ,  over  fairly  wide  temperature  ranges,  under  adiabatic  conditions, 
with  values  of  y  between  1.3  and  2  (see  [4]). 

For  small  deformations,  one  can  fit  (2)  adequately  by  the  lineari¬ 
zation  (y  =  -  1  ) 

P  2 

(4')  P  =  PQ  -  b  o,  F'(cr)  =  b  >  0,  a  constant. 

This  approximation  (Hooke's  Law)  is  used  in  linear  elasticity  theory; 
b  =  pc  =  c/a  is  the  rate  at  which  material  is  passed  by  a  sound  wave, 
c  being  the  sound  speed.  In  fluid  mechanics,  it  is  called  the  Chaplygin- 
Karman-Tsien  approximation.  Note  that  it  reduces  ( 1 )  to 


the  ordinary  wave  equation.  Alone  among  equations  of  state,  (4')  permits 
plane  waves  of  finite  amplitude  to  travel  without  change  of  form.  Cor¬ 
respondingly,  it  avoids  the  formation  of  shocks  [22,  §§  12,  1 33 •  Unfor¬ 
tunately,  it  is  not  satisfied  by  any  known  real  fluid. 
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In  classical  continuum  mechanics,  equations  (l)-(2)  are  assumed  to 

be  applicable  also  to  "slab"  motions  of  solids;  we  shall  discuss  this 

v  assumption  in  Part  II. 

Radial  cylindrical  and  spherical  motion  leads  to  analogous  equa- 
**  * 

tions.  By  using  appropriate  conditions  of  interfacial  continuity  in  ve¬ 
locity  u  =  x,  and  pressure  p,  one  can  treat  problems  involving  several 
homogeneous  layers  (slabs)  or  shells  in  contact  with  each  other.  And  by 
letting  F(a)  =  F(a,cr),  one  can  also  treat  problems  of  non-homentropic 
isentropic  flow. 

3.  Spatial  discretization.  In  order  to  solve  initial  value  (Cauchy) 

problems  for  (l)-(2)  on  a  digital  computing  machine,  one  must  use  a  finite 

set  of  values  aQ, ...,8^  of  the  Lagrangian  mass-variable  a,  and  a  finite 

set  of  times  t  <  t,  <  t_  <  •••.  One  can  then  think  of  x(a.,t.)  as  re- 
o  1  2  1  j 

presenting  the  position  of  the  i-th  particle  at  time  t ..  In  this  way, 

J 

each  finite  difference  approximation  to  ( 1 ) — ( 2)  can  be  regarded  as  defin¬ 
ing  a  law  of  motion  for  a  discrete  system  of  particles,  movements  being 
in  discrete  jumps. 

It  is  very  suggestive  to  fix  the  number  of  particles,  and  to  let 

At.  =  t.  -  t.  ,  tend  to  zero.  In  the  limit,  one  obtains  a  system  of  or- 
J  J  J-1  — 

dirmry  DE's  from  the  given  system  of  partial  DE's,  of  exactly  the  type 
(3)  considered  in  molecular  models  of  matter.  This  gives  a  spatial 
discretization,  sometimes  also  called  semi -discrete.  Such  spatial  dis¬ 
cretizations  have  been  considered  since  the  time  of  Lagrange  by  many 

*[5,  §  18];  the  first  exponent  2  in  the  last  formula  there  should  be 
deleted. 
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authors,  who  have  had  four  distinct  purposes  in  mind. 

First,  the  system  (3)  and  its  generalizations  have  been  considered 
as  defining  molecular  models  of  solid  crystals  by  Cauchy,  Kelvin,  Born, 
L.  Brillouin,  and  others  ([1],  [2],  [3h  [2*0).  These  authors  have  gen¬ 
erally  considered  the  small  deformations  typical  of  the  theory  of 
elasticity;  for  such  small  deformations,  linear  stress-strain  relations 
like  (41)  are  adequate. 

Second,  three-dimensional  analogs  of  (2)  have  been  applied  to  the 
kinetic  theory  of  dilute  gases,  by  Maxwell,  Boltzmann,  and  their  suc¬ 
cessors  ([4],  [8],  [17]).  The  force-laws  assumed  here  have  been  highly 
non-linear,  the  forces  being  negligible  except  occasionally,  during 
binary  encounters  corresponding  to  near-collisions.  These  encounters 
are  then  studied  statistically. 

Third,  analogs  of  (3)  have  been  applied  to  simulate  equations  of 
state,  and  especially  the  change  of  phase  from  dense  gases  to  liquids. 

A  resume  of  the  status  of  such  applications  may  be  found  in  [17,  Ch.  4]* 

Finally,  there  is  the  tradition  stemming  from  von  Neumann  [9],  in 
which  the  numerical  integration  of  (3)  is  used  to  compute  approximate 
solutions  to  the  non-linear  wave  equation  (l)-(2).  Since  a  character¬ 
istic  feature  of  this  work  consists  in  the  use  of  non-linear  force-laws 
which  are  almost  never  negligible,  it  corresponds  mathematically  most 
closely  to  molecular  models  of  dense  gases  and  liquids. 

*Such  models  were  first  studied  by  Newton,  who  speculated  that  the  phe¬ 
nomena  of  Nature  •  • .  may  all  depend  on  certain  forces  by  which  the 
particles  of  bodies  ...  are  either  mutually  impelled  towards # one 
another,  and  cohere  in  regular  figures,  or  are  repelled  ^lncjjja, 
preface  to  the  first  edition).  See  also  ibid.,  Book  II,  Prop. 
where  "near  neighbor"  force-laws  are  considered. 
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4.  Shocks;  equation  of  state.  It  is  well-known  that  the  homentro- 
pic  equation  of  state  (2)  breaks  down  physically  across  (strong)  shock 
fronts  in  gases  and  liquids.  When  such  shocks  occur,  one  must  replace 
(2)  by  a  thermodynamic  equation  of  state  of  the  form 

(5)  p  =  p(cr,T),  T  =  temperature. 

The  essential  feature  of  (5)  is  its  use  of  two  variables  in  place  of  one 
to  specify  the  local  state  of  matter. 

When  heat  conduction  k  and  viscosity  b  are  negligible  (adiabatic 
flow),  so  that  the  entropy  S  is  constant,  one  can  deduce  (2)  by  combining 
appropriate  assumptions  about  the  internal  energy 

(5*)  I=l(cr,T) 

with  energy  conservation  principles.  Thus,  for  a  perfect  gas,  I  =  c  T 
by  definition  and  (4)  follows  with  7  =  c  /c^.;  see 

By  combining  the  preceding  considerations  with  mass  and  momentum 
conservation,  and  allowing  heat  conduction  across  shocks,  one  can  then 
deduce  the  Rankine-Hugoniot  shock  relations  [5,  §§  5^-62].  These  are 
independent  of  the  specific  conductivity  k  assumed,  which  merely  affects 
the  thickness  of  the  shock  layer. 

It  is  generally  believed  that  the  inviscid  compressible  flow  equa¬ 
tions  of  1  1 ,  when  taken  with  the  Rankine-Hugoniot  equations  across 
shocks,  give  rise  to  a  system  of  partial  differential  equations  and  inter¬ 
face  conditions  which  define  a  well-set  initial  value  problem.  Moreover,  it 
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is  believed  that  the  resulting  system,  if  integrated,  adequately  simu¬ 
lates  real  fluid  motions  under  many  conditions.  However,  this  has  never 
been  rigorously  proved. 

It  is  further  believed  that  the  preceding  model  is  adequate  for 

* 

the  treatment  of  strong  shocks  in  solids.  If  supplemented  by  chemical 
kinetics  considerations,  it  can  be  used  to  treat  inviscid  hypersonic 
fluid  motions.  And  if  supplemented  by  realistic  fluid  friction  and  heat 
transfer  relations,  it  can  be  used  to  treat  fluid  flow  in  a  heated 
channel. 

5.  Viscous  compressible  flow^  shock  thickness.  However,  the  pre¬ 
ceding  hybrid  model  is  not  logically  satisfactory.  It  is  more  logical, 
and  more  consistent  with  physical  reality  generally,  to  treat  explicitly 
the  effects  of  viscosity  p  and  thermal  conductivity  k.  It  is  generally 
believed  that  the  model  of  §  4  can  be  obtained  rigorously  from  the  re¬ 
sulting  system  of  partial  differential  equations  as  the  limit  p  4,0, 
k  1 0. 

This  being  assumed,  it  seems  likely  that  the  limit  p  j,  0  with  k  =  0 
will  do  as  well  and  be  simpler.  It  is  primarily  this  limit  which  will 
be  discussed  below.  Only  for  a  quantitative  analysis  of  shock  wave 
thickness  is  k  important. 

We  therefore  treat  now  the  one-dimensional  viscous  compressible 

* 

See  [19];  shear  strength  is  perhaps  negligible. 

See  [5,(63.04)],  where  however  bulk  viscosity  is  neglected. 


flow  of  a  perfect  gas  in  the  limiting  case  k  =  0;  this  is  the  only  case 
for  which  computations  are  available. 

In  a  perfect  gas,  the  entropy  S  =  cy  In  (p  o7).  One  can  define  the 

"viscous"  stress  q  =  4  ^  \xj3  [18,(5)],  where  |i  =  |i(T)  is  a  material  con- 
* 

stant  depending  on  the  temperature.  The  one-dimensional  Euler- Lagrange 
equation  of  motion  (l)  is  then  replaced  by  the  one-dimensional  Navier- 
Stokes  equation 

(6)  xtt  =  ut  =  -  (P  -  <l)a,  q  =  m  ux  =  ua  =  HXta/xa. 

At  the  same  time,  the  polytropic  equation  of  state  (4)  is  replaced  by 
the  thermodynamic  equation  of  state  of  a  perfect  gas  [18,(4)] 

S/c 

(7)  P  =  c  (y  -  1 )T/o  =  e  v  d"7,  a  =  x  . 

V  Q, 

The  rate  of  change  of  entropy  is  [18,(5)] 

(8)  St  =  qxat/T  =  qu/T. 

Shock  thickness.  In  a  viscous  compressible  fluid,  shock  waves 
have  a  finite  thickness.  The  shock  thickness  can  be  estimated  by  assum¬ 
ing  time-independent  flow  in  (6)-(8).  Given  the  upstream  state  (S^,a^) 

** 

and  associated  sound  velocity,  there  exists  for  each  u^  >  c^  a  time- 
independent  flow  (unique  up  to  translation)  for  which  lim  u(x)  =  u, 

*Both  [18]  and  [5,§  63]  ignore  the  difficult  problem  of  treating  "bulk 
viscosity"  correctly;  for  this  problem,  see  [22,  §  33] • 

H.  Weyl,  Comm,  pure  appl.  math.  2  (1949),  103-22;  D.  Gilbarg,  Am.  J. 
Math.  73  (1951),  256-74. 


-15- 


lim  a(x)  =  a, .  For  this  flow,  u^  =  lim  u(x),  cr  =  lim  o(xl 

and  Sg  =  lim  S(x)  exist  and  are  given  by  the  Rankine-Hugoniot  equa¬ 
tions,  For  shocks  of  finite  amplitude  (with  cr^/o^  >  "*  *  say)>  most  of 

the  flow  is  nearly  uniform,  the  "shock  thickness"  6  over  which  99$  of  the 
change  in  cr(x)  takes  place  being  very  small.  The  shock  thickness  6  de¬ 
creases  with  the  strength,  being  of  the  order  of  the  molecular  mean  free 
path*  if  G^/tTg  >  2. 

6.  Real  and  artificial  viscosity.  If  a  physically  real  viscosity 
is  used  in  Eqs.  (6)-(8),  the  effect  is  negligible  except  in  shock  layers. 
Moreover  these  are  so  thin  in  comparison  with  mesh- lengths  of  interest 
that  real  viscosity  effects  cannot  be  computed  effectively  on  digital 
computing  machines,  in  typical  high-speed  flows  (large  Reynold  numbers). 

For  somewhat  the  same  reason,  it  is  also  awkward  to  treat  the  model 
of  i  4  on  computing  machines  with  a  rectangular  space-time  mesh.  The 
shocks  go  between  mesh-points  and  "shock-fitting"  must  be  done  by  inverse 
interpolation.**  To  indicate  what  is  done  in  practice,  we  write  down  the 
spatial  discretization  of  the  partial  differential  equations  (6)-(8). 

This  is  (since  pt  =  ([7  -  1]  q  -  7P>  these  equations): 

(9)  «n  -  -  (pn+|  -  Pn_i)  +  (qn+i  -  m  =  an  ‘  an-1' 

(10)  #a+A  -  ((7  -  1)^4  -  »n+J>  (*n+1  -  y/tVl  -  **)' 

*D.  Gilbarg  and  A.  Paolucci,  J.  nat.  mech.  anal.  2  (1953) 7  617-42.  See 
also  R.  Becker,  Z.  Phys.  8  (1922),  321-62;  L.  H.  Thomas,  J.  Chem.  Phys. 

12  (1944),  449-53;  H.  Grad>  Comm,  pure  appl.  math.  5(1952),  257-300. 

** 

See  L.  H.  Thomas,  Comm,  pure  appl.  math.  7  (1954), 195-206. 
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(11) 


a  i  =  —  p(3c  -  &  )/(x  -  -  x  ). 

^n+5  3  '  n+1  n'' v  n+1  n' 


The  pressure  relation  corresponds  to  the  equation  of  state  (7): 


(12) 


S/e 


P  1  =  e 


t(Vi  - 


Since  (9)-(l2)  define  a  system  of  ordinary  differential  equations,  it  is 
easy  to  integrate  the  system  numerically  with  high  accuracy  by  any  of  var¬ 
ious  methods. 

The  spatial  discretization  (9)-(l2)  reduces,  in  the  limit  p  =  0 
(q  =  0),  to  von  Neumann's  mass-spring  system  (3).  The  viscous  or  damping 
term  q  is  equivalent  to  inserting  a  dashpot  in  parallel  with  each  spring. 
This  model  was  proposed  for  solids  by  Kelvin  and  Voigt;  it  is  the  simplest 
model  of  rheology. 

To  achieve  smooth  velocity  profiles  with  a  finite  mesh- length,  one 
must  introduce  an  artificial  viscosity  large  enough  to  make  the  thickness 
of  the  shock  layer  two  or  three  mesh-lengths,  but  not  much  larger  than 
this.  To  achieve  this,  one  postulates  that  q  in  (6)  and  (11)  shall  de¬ 
pend  on  the  velocity  gradient  according  to  some  arbitrary  law  q(u  ,T) 

A 

(see  [6,  Part  IX],  [11],  [18],  [20]);  it  may  even  be  non-linear  in  ux. 

7.  Molecular  and  continuum  models  of  matter.  The  idea  that  con¬ 
tinuum  models  of  matter  can  be  derived  from  molecular  models,  as  the 
limiting  case  N  -»«  of  N  particles  of  equal  mass  m  =  Aa=  afc+1  - 

See  [12,  Ch.  X,  1  9]j  R.  Landshoff,  Report  LA-1950*  Lax  [20]  advocates 
using  a  linear  viscosity.  One  can  postulate  a  wave-front,  and  then  de¬ 
duce  the  associated  viscosity  "law"  (inverse  method). 
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very  old.  It  can  be  justified  for  elastic  fluids  as  follows 


If  the  molecules  attract  each  other  according  to  a  force-law 
F([xk+1  -  x^J/m),  the  equations  of  motion  are 


(13) 


mx 


1  =  F((xi+1  “  xj_Vm)  “  1  =  1/2 ,•••, 


N. 


Let  y  =  y(a,t)  be  a  function  with  two  continuous  derivatives  with  re¬ 
spect  to  a  such  that  y(a..,t)  =  x±(t),  i  =  1,2,***,N.  Then,  expanding  F 
in  Taylor  series  about  dy/da  for  a  =  a^,  substituting  into  ( 1 3)>  and  di¬ 
viding  by  m,  we  obtain 

(lit-)  dyta^tj/dt  =  F*(^y[ai,t]/2a)82y(ajL,t)/Sa2  +  o(Aa). 

In  the  limit  Aa  =  m  -»0  (N  -><»),  tends  to  (l)  if  ye  €  varies  sufficiently 
smoothly  in  position  and  time,  and  Fe  #  . 

The  mass-spring  model  ( 1 3)  and  its  generalizations  have  been  con¬ 
sidered  for  many  hypothetical  force-laws,  and  it  is  interesting  to  con¬ 
sider  some  peculiarities  of  these  laws.  We  begin  with  the  force-law 
f(r)  =  kr^  corresponding  to  the  polytropic  equation  of  state  (4). 

In  this  case,  the  particle  model  defined  by  ( 1 0)  admits  the  energy 
integral,  constant  in  time, 

05)  E  =  ^  m  £  x^  +  ^SL_  £  K*,  “  x^_i  7/1. 

(A  similar  integral  holds  for  any  F(d).) 

If  7  >  1  (a  logarithmic  analog  of  (15)  exists  for  7=1),  then  it 
takes  infinite  energy  to  make  two  particles  collide.  As  a  corollary,  the 
order  of  particles  cannot  change:  each  particle  always  has  the  same  two 
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"nearest  neighbors".  On  the  other  hand,  if  y  <  1 ,  it  only  takes  finite 
energy  to  make  x^  =  x_.  ^ ,  and  particles  can  collide  or  pass  "through" 
each  other.  (For  y  =  -1  and  for  particles  of  equal  mass,  such  inter¬ 
penetration  is  mathematically  equivalent  to  elastic  collision. ) 

It  is  curious  that  the  condition  y  >  1  should  also  coincide  with 
the  condition  that  c2  =  dp/dp  =  ky  should  increase  with  the  den¬ 

sity  •—  i.e.,  that  the  fluid  should  be  normal.*  (H.  Weyl  (op.  cit. 

2  2 

supra,  p.  105)  postulates  the  related  but  weaker  condition  that  p  c 
should  increase  with  the  density.)  Note  that  k/  >  0  since  F'(cr)  >  0  in 
(2). 

8.  Wave  transmission.  A  fundamental  description  of  real  and 
synthetic  homogeneous  materials  concerns  their  mode  of  transmission  of 
plane  waves  of  given  frequency  and  amplitude.  For  acoustic  or  small- 
amplitude  waves,  all  elastic  fluid  continua  are  characterized  by  their 
sound  velocity  cq  =  b/pQ;  the  general  wave  form  is  x  =  f(a  -  bt) 

+  g(a  +  bt)  in  Lagrangian  coordinates,  where  b  =  F'(oq).  In  Eulerian 
(space-time)  coordinates,  the  general  wave  form  6x  =  f(x  -  cQt) 

+  g(x  +  cQt)  is  the  same,  because  of  the  small  amplitude  approximation. 

Dispersion.  Using  Fourier  analysis,  one  can  attribute  the  preced¬ 
ing  result  to  the  fact  that  sinusoidal  acoustic  waves  travel  in  elastic 
fluid  continua  without  attentuation  or  dispersion.  This  is  however  not 


*See  [5,p.  5];  also  the  footnote  on  p.  5  of  G.  Birkhoff  and  J.  M.  Walsh, 
Riabouchinsky  Jubilee  Volume,  Paris,  1954.  Glass  may  be  "abnormal": 
dp/p^dp  seems  to  increase  with  p  up  to  10  kilobars. 
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true  in  molecular  models:  though  there  is  no  attenuation  (in  pass-hands), 
there  is  always  dispersion.  The  wave  velocity  c  of  a  sinusoidal  acous¬ 
tic  wave  of  length  A,  in  a  one-dimensional  mass-spring  system  with  par¬ 
ticles  spaced  h  =  l/N  apart,  satisfies  [3,  p.  4] 

(16)  c  =  c  sin  ( jth/A)/(jth/A), 

00 

where  c  is  the  velocity  of  long  waves  (A  -*<»).  That  is,  there  is 
00 

dispersion:  c  varies  with  A,  between  the  limits  c^  and  2c Jx  ~  0.63350^. 

Though  the  preceding  phenomenon  of  dispersion  is  not  observed 
physically  in  gases,  it  does  simulate  (with  molecular  particle  spacing) 
some  characteristic  features  of  the  behavior  of  solid  crystals.  For 
example,  many  crystals  have  a  frequency  cutoff..  Other  models  consisting 
of  two  kinds  of  masses,  alternatively  spaced,  will  transmit  waves  only 
in  certain  frequency  bands.  In  this  respect,  such  models  simulate  the 
optical  behavior  of  many  crystals  and  the  behavior  of  electric  band-pass 
networks;  see  (3]. 

For  large-amplitude  waves  in  elastic  fluid  continua,  the  wave  form 
is  still  given  in  Lagrangian  coordinates  by  x  =  f(a  -  cQt)  +  g(a  +  CQt) 
in  the  case  y  -  -1  of  linear  elasticity  theory  (Hooke's  Law),  with  which 
most  of  the  relevant  literature  has  been  concerned.  In  molecular  models, 
(16)  still  holds. 

Shocks.  It  is  well-known  that  shocks  arise  in  polytropic  elastic 
fluids  if  y  >  -1 .  Moreover  simple  shocks  (with  uniform  flow  on  each 
side  of  the  shock)  are  transmitted  without  change  of  form  for  any 
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(normal)  equation  of  state  in  any  continuous  elastic  fluid.  (For  simple 
shocks  in  viscoelastic  media,  see  §5.) 

As  a  result  of  dispersion,  this  feature  is  not  reproduced  in  the 
spatial  discretization  discussed  above,  even  when  y  =  -1 .  This  was 
first  observed  by  von  Neumann  [9],  for  y  =  1 .4  (approximately). 

When  y  =  - 1,  the  analog  for  molecular  models  of  Kelvin's  Principle 
* 

of  Stationary  Phase  is  applicable  to  simple  shocks  and  other  pulses  of 
short  duration.  This  principle  asserts  that  each  Fourier  component 
travels  with  its  own  wave  velocity,  so  that  the  appearance  of  a  pulse  of 
short  duration  after  being  transmitted  through  a  long  distance  should  be 
as  follows.  The  "head"  of  the  pulse  (which  travels  most  rapidly)  should 
largely  conserve  its  initial  form;  but  the  "tail"  should  be  wavy  with 
wavelength  diminishing  continuously  to  the  shortest  observable  wave¬ 
length  2/N.  It  would  be  interesting  to  test,  by  numerical  calculation, 
whether  this  was  also  true  in  the  non-linear  case,  with  y  >  -1 . 

One  can  speculate  that  dispersion  will  also  make  shock  wave  thick¬ 
ness  in  molecular  models  greater  than  in  continuum  models  (see  I  5) •  It 
would  be  interesting  to  study  this  question  for  the  viscoelastic  molecu¬ 
lar  model  of  1  6,  looking  for  solutions  periodic  in  time.  If  the  equa¬ 
tions  were  linear  (7  =  -l),  the  shock  thickness  would  increase  indefi¬ 
nitely  with  time. 


*[3,  p.77];  [23,  §  2kl];  J.  J.  Stoker,  "Water  waves",  Interscience  (1957)> 
p.  163. 
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9.  Finite  time-steps.  In  §  1-8,  we  have  discussed  systems  of 
particles  whose  positions  were  continuously  defined  in  time,  of  a  type 
familiar  in  the  kinetic  theory  of  gases  and  in  the  molecular  theory  of 
crystals.  In  digital  computations  one  uses  finite  time-steps  At. 

It  is  important  to  choose  these  At  so  as  to  minimize  the  computa¬ 
tional  (truncation  and  roundoff)  errors.  In  the  special  case  y  -  -1, 
both  truncation  and  dispersion  can  be  eliminated  by  using  At  =  Ax/c, 
i.e.,  by  integrating  along  characteristics .  The  accuracy  of  the  approxi¬ 
mation  depends  solely  upon  roundoff  and  the  accuracy  of  determining  x^O) 
and  x^(At  ). 

Convergence  theorems.  The  discussion  at  the  beginning  of  §  7  sug¬ 
gests  that,  as  Aai  tends  to  zero,  the  solutions  of  suitable  difference 
approximations  to  (l)— (2)  should  converge  to  the  solution  of  the  differ¬ 
ential  equations,  for  the  same  initial  conditions.  The  natural  conjec¬ 
ture  is  that  this  is  true  provided  no  shock  waves  occur  and  that  the  time- 
steps  Ati  satisfy  the  stability  condition  At±  <  Ax/c,  the  local  sound 
speed. 

For  the  differential  equations  describing  the  motion  of  an  elastic 
fluid  in  Eulerian  coordinates, 

ut  +  uux  =  -  c  (p)  p  Px,  Pt  +  (pu)x  *  0, 

convergence  follows  from  a  general  theorem  of  Courant,  Isaacson  and  Rees. 

It  also  follows  for  the  coupled  first-order  hyperbolic  equations 

*R.  Courant,  E.  Isaacson,  and  M.  Rees,  Comm,  pure  appl.  math.  2  0952), 
24>55. 
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ut  =  b2(a)ca,  ua  =  at,  b (a)  =  Vf'(ct), 

equivalent  to  (l)-(2),  from  the  same  general  theorem.  To  obtain  x(a,t) 
from  the  approximate  solution,  one  must  then  integrate  the  equation 
xa  =  a. 

The  natural  difference  approximation  to  the  second-order  hyperbolic 
equation  (l),  is 

(17)  Xj  -  2Xj  +  x^.  =  “Si,  “  *j]/Aa)  ~  F([xj  “  xj_iVA^}> 

where 

Xj  =  x(jAa,  nAt), 

* 

The  solution  of  this  set  of  difference  equations  converges  to  the  solu¬ 
tion  of  (l)  as  At  and  Aa  go  to  zero,  subject  to  the  stability  requirement 
At^F ' /Aa2  <  1;  this  requirement  is  met  if  At/Aa  is  chosen  small  enough  so 
that  Ax^/At  is  greater  than  the  local  speed  of  sound 

”  _i  i_ 

c  =  a|dp/da|2  =  o|f'(c)|2. 

As  a  corollary,  letting  first  At^O  for  fixed  Aa  and  then  letting 
Aa-1'0,  one  obtains  a  rigorous  proof  for  the  convergence  of  the  spatial 
discretization  ( 1 )— (2)  to  the  exact  solution  of  the  corresponding  con¬ 
tinuum  problem,  provided  there  are  no  shocks. 


2  2  1 

Provided  F c-6  ,  x(a,0)e#  ,  and  x  (a,0)e-#  ,  and  that  the  derivatives 
Fff(J,x  (a,0)  and  xat(a,0)  satisfy  Lips  chit  z*  conditions.  The  proof  will 
be  given  in  the  thesis  of  Robert  E.  Lynch. 
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Unfortunately;  no  similar  results  have  been  proved  for  the  one¬ 
dimensional  Navier-Stokes  equations  (6)-(8).  Even  the  convergence  of 
solutions  of  (9)-(H)  to  solutions  of  (6)-(8)  as  m  -*  0  (N  -+«)  has  not 
been  established.  Neither  has  any  convergence  theorem  been  proved  in 
the  limit  Aa  -* 0,  p  -» 0  to  the  motion  of  an  elastic  fluid  with  shocks. 

Synthetic  materials.  Referring  back  to  §  1,  it  is  interesting  to 
speculate  as  to  what  kind  of  synthetic  materials  are  defined  by  computa¬ 
tion  schemes  like  ( 1 7)  with  finite  time-steps  At  .  This  question  can  be 
asked  also  for  difference  approximations  to  the  equations  of  a  visco¬ 
elastic  fluid.  For  some  computing  processes  (such  as  forward-difference 
methods),  the  use  of  discrete  time-steps  can  be  regarded  as  introducing 
a  finite  delay  time  into  reactions.  Such  delay  times  or  relaxation  times 
have  also  been  suggested  as  being  responsible  for  various  physical  phe- 
nomena,  such  as  ultrasonic  attenuation.  In  the  theory  of  viscoelastic 
solids,  similar  delay  times  characterize  Maxwell  solids.  Their  charac¬ 
teristic  property  is  to  cause  high  attenuation  in  narrow  resonance  bands 
—  as  contrasted  with  monotonically  varying  viscous  attenuation.  Such 
resonance  bands  are  observed  physically. 

Further  study  may  reveal  other  qualitative  phenomena  associated 
with  the  use  of  finite'  time- steps,  also  observable  physically.  By  study¬ 
ing  such  phenomena  in  one  space  dimension,  one  can  hope  to  gain  insight 
into  the  three-dimensional  case. 

*See,  for  example,  M.  Brillouin,  "La  viscosite  des  liquides  et  des  gaz," 
Paris,  1907,  vol.  2,  pp.  93-5;  K.  F.  Herzfeld  in  Part  H  of  "Thermodyna¬ 
mics  and  physics  of  matter,"  F.  D .  Rossini,  editor,  Princeton,  1955* 
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II.  STATISTICAL  MECHANICS 


10.  Kinetic  theory  of  gases.  In  the  kinetic  theory  of  gases,  it 

* 

is  assumed  that  gases  are  composed  of  very  many  moving  particles,  inter¬ 
acting  with  each  other  according  to  appropriate  molecular  force-laws.  In 
the  semi-discrete  model  (§  3)  of  a  shock  advancing  into  a  uniform  gas,  ir¬ 
regular  particle  motion  occurs  behind  the  shock;  von  Neumann  suggested 
[9]  that  this  motion  was  analogous  to  the  molecular  motion  in  the  kinetic 
theory  of  gases.  We  will  now  analyze  this  suggestive  idea  critically. 

In  the  kinetic  theory  of  gases,  pressure  p  and  temperature  T  are 
independent  variables,  as  in  the  continuum  model  of  §  4.  Moreover,  the 
temperature  has  a  purely  mechanical  significance,  as  the  mean  transla¬ 
tional  energy  of  molecular  motion.  Thus,  by  analogy  with  the  three- 

** 

dimensional  case,  we  can  write 

(18)  RT  =  1  m(x  -  x)2  =  ^  Z  mi(xi  -  |  Z  x^2, 

allowing  for  the  possibility  of  variable  mass  m^.  Here  x^  is  the 

About  N  =  2.7  x  10  y  per  cc  (Avogadro’s  number)  under  standard  atmos- 

07 

pheric  conditions  [8,  p.  8];  the  incorrect  value  6x10^  (Loschmidt's 
number) is  frequently  quoted. 

See  [8,  §  146,  p.  31l],and  compare  with  §  l4l,  p.  299  and  §  26,  p.  32. 
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position  of  the  i-th  mass  point  in  the  semi -discrete  model  of  §  15  the 
summations  are  over  N  gas  molecules.  (Note  that,  by  Maxwell's  Principle 
of  the  Equipartition  of  Energy  [8,  p.  8l],  the  statistical  expectation 
of  E^  =  nu (x_^  -  x)  /2  is  independent  of  i.)  Whereas  in  mass-spring 
models  of  elastic  fluids,  which  reproduce  the  adiabatic  equation  of  state 
under  quasi-static  conditions  (i.e.,  very  gradual  compression),  the  pres¬ 
sure  p  determines  the  nominal  temperature  T. 

As  N  -*  os,  the  convergence  theorems  of  I  9  show  that  motions  of 
most  individual  molecules  in  the  mass-spring  model  (3)  will  deviate  less 
and  less  from  the  mean  motion,  except  behind  shocks.  That  is,  the  effec¬ 
tive  local  temperature  (in  the  sense  of  kinetic  theory)  is  zero.  In  par¬ 
ticular,  the  intermolecular  forces  f(r)  between  adjacent  particles  in 
Lagrangian  hydrodynamical  computations  act  continuously  and  normally  vary 
by  a  factor  of  at  most  ten,  whereas  the  intermolecular  forces  treated  in 
the  kinetic  theory  of  dilute  gases  are  negligible  most  of  the  time,  being 
appreciable  only  during  brief  near- collisions  (encounters). 

Specific  heat.  The  proportionality  of  temperature  to  the  kinetic 
energy  of  relative  motion  (relative  to  the  local  mean  motion)  is  clearly 
illustrated  by  a  study  of  specific  heats.  In  monatomic  gases,  the  speci¬ 
fic  heat  is  about  3  cal/mole  °C  [4,  p.  42]  above  the  critical  point;  in 
diatomic  gases,  it  is  about  5  cal/mole  °C.  In  solids,  where  the  poten¬ 
tial  energy  of  vibration  and  the  kinetic  energy  of  vibration  are  about 
equal,  it  is  often*  about  6  cal/mole  °C.  The  number  3  is  equal  to  the 

*This  is  the  Dulong  and  Petit  Law  [3,  p.  1 66] ,  approximately  valid  if  the 
"Debye  temperature"  is  one  or  greater. 
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number  of  degrees  of  freedom  of  motion,  due  to  a  lucky  coincidence: 
that  H^O,  with  molecular  weight  18,  has  the  abnormally  large  specific 
heat  of  18  cal/mole  °C  at  20°C. 

At  low  temperatures,  the  specific  heats  of  solids  are  consid¬ 
erably  less  than  6  cal/mole  °C.  For  various  force-laws  f(r),  it  would 
be  interesting  to  study  the  mean  energy  in  thermal  equilibrium  of  a 
"solid"  consisting  of  n  particles  of  mass  m,  bombarded  at  random  times 
by  particles  of  mass  m  having  a  Maxwell  velocity  distribution  and  mean 
energy  mv  /2.  From  this  one  could  infer  the  ratio  of  specific  heats 
and  the  Debye  temperature  curve  for  various  one-dimensional  "solids". 

11.  Equation  of  state.  In  the  mass-spring  model  of  §  j>,  the 
pressure  is  attributed  to  static  molecular  repulsion;  in  the  kinetic 
theory  of  gases,  it  is  attributed  to  molecular  motion.  The  first  model 
simulates  liquids  and  solids,  to  a  first  rough  approximation,  while  the 
second  simulates  dilute  gases  very  well,  especially  monatomic  gases. 

That  is,  pressure  in  gases  is  dominantly  kinetic,  whereas  in  liquids 
and  solids  it  is  dominantly  static. 

In  general,  the  virial  of  Clausius  allows  one  to  decompose  the 
pressure  into  a  kinetic  component  proportional  to  RNT,  and  a  static  com¬ 
ponent  proportional  to  Nrf(r) ,  where  N  is  the  number  of  particles  per 
unit  volume.  In  n  dimensions, 

(19)  per  =  ^  [RET  +  Nrf[r)]  =  ^  [RT  +  rf(r)]. 

See  [8,  §§  1 63-733 y  where  the  three-dimensional  case  is  treated. 
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Here  T  is  defined  as  the  mean  kinetic  energy  of  translation  per  particle. 

One  of  the  oldest  problems  in  the  molecular  theory  of  matter  is  to 
simulate  the  equations  of  state  of  real  material  by  molecular  models, 
postulating  an  appropriate  intermolecular  force-law.  Among  the  force- 
laws  which  have  been  studied  ([4],  [8],  [17]  ),  ^he  hardsphere,  power  laws 
f(r)  =  Ar"a,  and  the  Lennard-Jones  formula  f(r)  =  Ar”a  -  Br  ^  have  re¬ 
ceived  the  most  attention.  The  problem  is  to  deduce  the  observed  pres¬ 
sure  p(c,T),  internal  energy  E(cr,T),  and  phase  transitions  as  functions 
of  the  specific  volume  a  and  temperature  T. 

Because  such  studies  ignore  the  fact  that  real  matter  is  made  up 
of  positively  charged  atomic  nuclei  and  negatively  charged  electrons  of 
much  smaller  mass,  their  conclusions  cannot  be  taken  as  having  a  direct 
physical  significance.  Rather,  they  must  be  considered  as  suggestive. 

In  particular,  they  may  suggest  computation  schemes  for  simulating  the 
behavior  of  real  materials  in  (say)  atomic  explosions  more  realistically 
than  would  be  possible  with  an  idealized  elastic  or  inelastic  fluid. 

For  such  purposes,  studies  in  one  space  dimension  have  much  to  re¬ 
commend  them.  Moreover  in  one  dimension,  it  is  easy  to  supplement  ana¬ 
lytical  studies  by  numerical  studies  of  the  motions  of  sample  systems  of 
particles.**  We  hope  that  such  studies  will  be  made.  Because  of 

*The  Lennard-Jones  formula  is  also  assumed  by  Griineisen,  in  an  attempt 
to  correlate  the  Debye  temperatures  of  solids  with  their  compressibili¬ 
ties;  see  Born  and  Huang  [2,  p.  52]. 

**Such  a  study  has  recently  been  made,  for  the  rigid-sphere  model  in 

three  dimensions,  by  Alder  and  Wainwright.  See  "Transport  processes  in 
statistical  mechanics,"  I.  Prigogine,  editor.  Interscience,  1958,  97-140 
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similarity  (see  §  14),  one  should  be  able  to  compute  some  quantities 

with  N  <  100.  In  [9],  von  Neumann  suggested  (in  another  context)  that 

N  =  l4n  might  be  adequate  in  n  dimensions. 

In  anticipation  of  such  studies,  we  make  a  few  elementary  remarks 

about  statistical  mechanics  in  one  space  dimension. 

In  the  rigid- sphere  model,  it  is  evident  that  the  equation  of  state 

2 

will  be  p(cr  -  a  )  =  NT,  where  T  =  me  /2  and  cr  =  Nr  is  the  one-dimensional 
*  '  o  7  '  oo 

"volume"  occupied  by  molecules. 

In  general,  one  must  distinguish  between  forces  assumed  to  act  only 
between  nearest  neighbors,  as  in  mass-spring  systems,  and  universal 
forces  (like  gravitation)  assumed  to  be  defined  for  all  particle-pairs. 

We  will  here  consider  only  forces  of  the  first  kind.  These  were  origi¬ 
nally  postulated  by  Newton  (loc.  cit.  in  I  3),  and  they  are  better  suited 
to  numerical  computation. 

In  general,  letting  f.  (t)  =  p,  i(t)  be  the  force  exerted  on  the  k-th 

K  K-  2 

particle  "by  the  (k-l)-st  particle,  it  is  easily  seen  that  fk+r 

averaged  over  a  long  period  of  time.  Otherwise,  the  k-th  particle  would 
undergo  a  net  acceleration.  Hence,  f  must  be  independent  of  k  in  the 
one-dimensional  case. 

Now  consider  the  case  y  =  -1  of  Hooke’s  Law.  In  this  case,  if 
f(r)  =  pQ  -  Br  is  the  repulsive  force,  then  f^  =  pQ  -  B(x^  -  x^_^). 

Since  f  =  f  is  independent  of  k,  it  follows  that  Nf  =  NpQ  -  B(x^  -  x^). 
Solving  for  f  =  pQ  -  ^(x^  -  x^)/N,  we  see  that  the  pressure  p  =  f  satis¬ 
fies 
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(20) 


P  =  PQ  -  Bm a,  o  «  -  Xi )/l»m. 

This  shows  that  the  pressure  is  independent  of  the  temperature  (kinetic 
energy).  In  (19),  the  term  in  rf(r)  is  decreased  by  molecular  motion 
just  enough  to  compensate  for  the  increase  in  RNT;  the  Clausius  decompo¬ 
sition  (19)  is  completely  misleading. 

Similar  considerations  suggest  that,  in  the  synthetic  gas  defined 
by  the  mass-spring  model  of  §  3>  the  total  pressure  will  increase  with 
the  temperature  (we  will  have  (dp/dT)^  >  0)  if  7  >  -1 ,  so  that  repulsion 
increases  more  than  linearly  with  approach  (case  of  "stiff"  springs). 

In  the  opposite  case,  7  <  -1  of  "soft"  springs,  (dp/dT)0  will  be  nega¬ 
tive  in  (3)»  The  simplest  model  of  elastic  collision  is  the  extreme 
case  of  a  "stiff"  spring;  water  below  4°C  corresponds  to  the  hypothesis 
of  a  "soft"  spring;  7  =  2  corresponds  to  an  inverse  square  repulsion 
law. 

It  would  be  desirable  to  study  the  dependence  of  p(o,T)  on  f(r) 
quantitatively  by  the  method  described  above.  Using  larger  computing 
machines  such  as  STRETCH  and  LARC,  one  could  try  to  simulate  more  gen¬ 
erally  the  statistical  mechanics  of  polar  molecules  and  the  elusive 

* 

phenomena  of  condensation  and  evaporation. 

Monte  Carlo  methods.  Monte  Carlo  methods  have  also  been  used  to 
calculate  the  equation  of  state  for  various  intermolecular  force-laws 

See  [17,  Ch.  5lj  H.  N.  Temperley,  "Changes  of  state,"  London,  1956, 
pp.  2-3;  0.  K.  Rice,  Part  E  of  "Thermodynamics  and  physics  of  matter," 
Princeton,  1 955 ^  C.  N.  Yang,  Phys.  Rev. 
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* 

—  and  especially  for  rigid  elastic  spherical  molecules.  This  approach 
is  related  to  that  discussed  in  this  report  only  in  its  objective  and 
its  use  of  high-speed  computing  machines. 

12.  Ergodic  theory.  Though  the  conditions  of  statistical  equili¬ 
brium  assumed  in  kinetic  theories  of  matter  are  not  attained  during 
times  typical  of  hydrodynamical  interest,  it  is  interesting  to  speculate 
what  these  conditions  are,  and  what  is  the  order  of  magnitude  of  the 
time  scale  required  to  attain  them.  High-speed  computing  machines  may 
be  expected  to  provide  useful  information  regarding  these  questions,  and 
hence  indirectly  regarding  molecular  models  of  matter  based  on  equations 
of  the  form  (3)  and  statistical  mechanics. 

In  using  computing  machines  to  perform  such  research,  it  is  impor¬ 
tant  to  realize  that  some  of  the  fundamental  concepts  of  statistical 
mechanics  have  themselves  never  been  rigorously  established.  For  in¬ 
stance,  it  has  never  been  proved  that,  as  t  »,  the  only  statistical  in¬ 
variants  of  an  assembly  of  molecules  in  a  rigid  container  with  elastic 
walls  are  temperature  and  density.  Stated  more  abstractly,  one  of  the 
gaps  in  the  theory  is  the  lack  of  proof  of  the  assumption  of  "metric 
transitivity"  (also  called  the  ergodic  hypothesis):  that  all  time- 
histories  having  given  energy  and  density  have  equivalent  statistical 
averages,  with  probability  one. 


The  pioneer'  study  (for  a  two-dimensional  "gas")  was  by  N.  Metropolis, 
A.  W.  and  M.  N.  Rosenbluth,  and  A.  H.  and  E.  Teller,  J.  chem.  phys.  21 
(1953),  1087-92.  For  subsequent  work,  see  W.  W.  Wood  et  al.,  Wuovo 
Cimento  9  0958),  1 33-^3;  Z.  W.  Salsburg  et  al.,  J.  chem.  phys.  30 
(1959),  67-72,  and  refs,  given  there. 
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Though  various  theoretical  arguments  tending  to  confirm  this  hypo- 
thesis  have  been  given,  it  still  rests  on  an  essentially  experimental 

basis,  and  the  generality  of  the  conditions  under  which  it  may  be  expected  * 

to  hold  is  not  clear.  Without  it,  special  arguments  are  needed  to  show 
that  the  function  p(cr,T)  is  well  defined  (i.e.,  single -valued) . 

And  yet  it  is  easy  to  show  that  the  ergodic  hypothesis  stated 
above  does  not  always  hold.  For  instance  the  power  spectrum  (energy  den¬ 
sity  as  a  function  of  wavelength)  is  an  invariant  in  a  Chaplygin  fluid 
with  y  =  -1 .  In  a  three-dimensional  cube,  the  energy  density  as  a  func¬ 
tion  of  the  wave-vector  ( |k^ |, Ik^J , |k^| )  is  invariant. 

Turning  to  one -dimensional  molecular  models,  it  is  classic  that 
for  N  rigid  elastic  spheres  in  one  dimension,  the  entire  set  of  veloci¬ 
ties  is  conserved.  Hence,  again,  density  and  temperature  are  not  a  com-  r 

plete  set  of  invariants,  even  though  p(o,T)  =  NET/ (a  -  <jq)  is  a  single- 

0 

valued  function. 

In  the  discrete  analog  of  a  Chaplygin  fluid,  defined  by  a  mass¬ 
spring  system  with  Hooke’s  Law,  there  is  no  energy  interchange  of  energy 
between  different  wavelength  (i.e.,  different  portions  of  the  spectrum). 

Hence  the  motion  is  not  metrically  transitive.  However,  since  the  fre¬ 
quencies  cn.  =  2  sin  (jit/2N)  of  the  normal  modes  are  ordinarily 
<3 


The  most  important  is  that  of  Oxtoby  and  Ulam,  Annals  of  Math.  kO  (1939) > 

560-66.  J.  Moser  has  recently  proved  a  complementary  unpublished  result 
on  the  existence  of  invariants.  For  a  physical  discussion,  see  G. 

Uhlenbeck  in  Appendix  I  of  Mark  Kac,  "Probability  in  the  physical  sciences," 

Interscience,  1959*  0 
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incommensurable*  (i.e.,  rationally  independent),  one  may  expect  phase 
randomization  by  dispersion  in  discrete  models. 

One  can  surmise  that,  in  the  non-linear  case  7  ^  -1 ,  both  fre¬ 
quency  and  phase  will  tend  to  distribute  themselves  statistically  in  a 
uniform  way,  so  that  time-averages  will  again  be  independent  of  initial 
conditions.  This  is  suggested  by  the  Theorem  of  Oxtoby  and  Ulam;  some 
interesting  experiments  bearing  on  it  have  been  performed  by  Fermi,  Pasta, 
and  Ulam  [1 6],  Some  further  experiments  along  similar  lines  for  a 
single  "ping-pong  ball"  have  been  reported  and  analyzed  at  the  Fourth 
Berkeley  Symposium  by  S.  Ulam  and  J.  M.  Hammersley. 

One  can  also  surmise  that,  in  two  or  more  dimensions,  molecular 
models  of  fluids  will  be  ergodic  for  arbitrary  force-laws  with  probability 
one  (i.e.,  for  almost  all  initial  velocities  and  distributions).  This  is 
because  of  the  great  variety  of  possible  angles  of  scattering.  (We  ex¬ 
clude  cases  like  that  of  closely  packed  rigid  spheres,  where  the  neigh¬ 
bors  of  a  fixed  "particle"  cannot  change.) 

13.  Diffusivities;  correlation  functions.  From  statistical  stu¬ 
dies  of  time-histories  of  assemblies  of  particles  subjected  to  various 
assumed  force-laws  f(r),  one  can  also  try  to  estimate  the  viscosity 


This  is  especially  likely  if  N  is  a  prime,  because  the  cyclotomic 
equation  is  then  irreducible  (Eisenstein's  Theorem). 

These  experiments  apparently  indicate  that  equipartition  of  energy 
does  not  occur.  Joseph  Ford,J.  Math.  Physics  2  (1961),  387-93*  has 
suggested  an  explanation  of  this.  His  suggestion  is  that,  for  the 
number  of  mass-points  used  in  the  experiments,  there  is  no  appreciable 
energy  sharing  between  the  weakly  coupled  oscillators. 
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(momentum  diffusion),  thermal  conductivity  k  (heat  diffusion),  and  ma¬ 
terial  diffusivity  D  associated  with  molecular  models  (j)  of  fluids.  If 
the  ergodic  hypothesis  is  fulfilled,  time-averages  should  he  independent 
of  the  initial  conditions. 

For  dilute  gases  (i.e.,  at  high  temperatures),  such  studies  will 
probably  not  be  immediately  rewarding.  Classical  methods  of  estimating 
binary  deflection  probabilities  seem  more  powerful  as  mathematical  tech¬ 
niques.  Moreover  physically,  the  ratio  of  the  three  diffusivities  (for 
momentum,  heat, and  matter)  seems  to  be  nearly  independent  of  the  force- 
law  assumed.  Thus  the  Prandtl  number  |i  C^/k  is  about  0.4  [4,  Ch.  10], 
while  Dp/p  is  between  1.0  and  1.5  [8,  Ch.  13]  for  spherically  symmetric 
molecules.  Deviations  of  real  gases  from  these  laws  seem  to  be  due 
largely  to  the  fact  that  real  gas  molecules  have  energy  of  rotation  as 
well  as  energy  of  translation,  unless  they  are  monatomic. 

However,  for  dense  gases  and  liquids,  one  can  probably  best  esti¬ 
mate  the  dependence  of  diffusivities  on  the  force-law  assumed  from  numeri¬ 
cal  experiments .  Such  experiments  might  also  reveal  interesting  rela¬ 
tions  between  the  shear  viscosity  p  and  the  bulk  viscosity  p',  which  is 
zero  in  a  dilute  monatomic  gas. 

To  get  the  most  information  out  of  such  (necessarily  elaborate) 
numerical  experiments,  a  great  deal  of  ingenuity  will  certainly  be  re¬ 
quired.  This  is  especially  true  of  dynamical  studies:  The  statistical 
mechanics  of  non-uniform  gases  is  very  difficult. 

For  example,  let  the  kinetic  temperature  T(x..)  of  a  mass-spring 


system  ("gas")  with  a  temperature  gradient  he  defined  in  statistical 
equilibrium  as  the  mean  kinetic  energy  nuu^/2  of  the  particle  with 
mean  position  =  x(a^) .  If  one  defines  the  kinetic  temperature  T(x,t) 
by  an  average  over  phase- space,  how  can  one  compute  this  effectively 
from  a  small  sample  of  a  few  nearby  particles?  What  is  the  probable 
error? 

To  obtain  some  information  on  this  question,  in  the  case  of  sta¬ 
tionary  "gases"  containing  just  one  kind  of  molecule  of  mass  m,  one  can 
define  a  closely  related  correlation  function 

(21)  0(k)  =  uiui+k  /  u\,  UjL  =  ij. 

This  will  provide  a  measure  of  the  smoothness  of  a  computation,  and  is 
also  a  fundamental  statistical  entity.  The  evolution  of  0(k,t)  with 
time  would  be  especially  interesting  to  study. 

The  study  of  such  correlation  functions  should  help  in  separating 
out  local  mean  motions  (those  which  one  would  obtain  in  the  limit  N  -»«») 
from  random  molecular  motions.  Hence,  indirectly,  it  would  help  in  the 
effective  definition  and  computation  of  conductivities  and  viscosities. 
By  analogy  with  the  theory  of  turbulence,  it  would  also  lead  to  the  in¬ 
teresting  concept  of  an  energy  spectrum  of  molecular  motion,  which  could 


*To  avoid  difficulties  associated  with  u  u.  k  being  undefined,  it  is 
inconvenient  to  use  a  periodic  gas  with  u^(t)  =  x^(t)  +  X. 

For  some  other  ambiguities  in  the  kinetic  theory  of  viscosity,  see 
C.  Truesdell,  Zeits.  fur  Physik  131  (1952),  273-89. 
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be  studied  as  a  function  of  the  force-laws  being  considered. 

14.  Springs  with  dashpots.  The  preceding  discussion  referred  to 
mass-spring  models,  with  emphasis  on  nearest  neighbor  force-laws  in  one 
space  dimension.  It  follows  from  the  discussion  of  I  13  that  such  mo¬ 
dels  will  have  a  molecular  viscosity  p.  Offhand,  one  would  expect  p  to 
be  of  the  order  of  magnitude  of  pvh/3,  where  v  is  the  mean  molecular  ve¬ 
locity  and  h  (the  particle  spacing)  is  the  analog  of  the  mean  free  path. 
It  would  follow  that  poo/T,  as  in  elementary  kinetic  theory. 

Of  great  interest  are,  also,  molecular  models  constructed  by  for¬ 
mulas  (9)-(l2)  from  the  Navier-Stokes  equations.  These  provide  molecu¬ 
lar  models  for  viscoelastic  fluids  and  for  viscoelastic  solids  of  various 
kinds.  These  models  also  describe  some  computation  schemes  using  an 
artificial  viscosity  p*  »  p  to  smooth  out  the  computed  flow  behind 
shocks.  In  such  computation  schemes,  where  the  pressure  depends  on  an 

internal  energy  (temperature)  T  for  given  a ,  one  can  introduce  an 
energy  conservation  law  (and  prove  that  the  entropy  increases  monotoni- 
cally).  The  temperatures  may  be  attributed  to  the  masses  or  the 
springs;  we  omit  the  formulas. 

*  [8,  Ch.  XI].  For  real  gases,  pccTn,  where  O.65  <  n  <  1. 

**(21,  Part  II] j  W.  P.  Mason,  "Physical  acoustics  and  properties  of  sol¬ 
ids,"  Van  Nostrand,  1958,  Ch.  VII;  D.  R.  Bland,  "Linear  viscoelasti¬ 
city,"  Pergamon  Press,  i960.  In  their  models  of  solids,  some  rheolo- 
gists  would  also  include  a  block  on  a  rough  surface  to  reproduce  the 
transition  to  plastic  flow  at  the  yield  point.  Springs,  dashpots, 
and  blocks,  in  various  series  and  parallel  combinations  (Bingham, 
Prandtl,  Kelvin,  Maxwell  bodies,  etc.)  have  been  found  useful  in  re¬ 
producing  characteristics  of  dough  and  other  substances. 
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For  each  choice  of  y,  artificial  viscosity  n*,  and.  spacing  h,  one 
can  thus  define  a  synthetic  gas,  whose  effective  viscosity  is  a  combina¬ 
tion  of  the  molecular  viscosity  and  the  artificial  viscosity  postulated. 
(This  is  analogous  to  combining  molecular  viscosity  with  turbulent  "eddy 
viscosity".) 

Similarity.  One  can  reduce  the  number  of  parameters  required  to 
define  such  synthetic  gases  by  identifying  systems  of  particles  whose 

behavior  is  similar  up  to  change  in  scale.  For  linear  viscosities,  one 

_  _  JL 

can  thus  define  a  natural  Reynolds  number  Re  =  p  ch/p*, where  c  oc  T2. 
Synthetic  "dilute"  gases  having  the  same  y  and  Re  should  be  completely 
similar. 

Analogous  remarks  apply  to  thermal  conductivity,  which  is  given 
quite  well  in  dilute  gases  by  the  formula  Pr  =  p  C^/k  =  4/(9 y  ~  5)  f°r 
Prandtl  number  [8,^  301]. 

Even  in  dense  gases,  the  addition  of  a  constant  pQ  to  the  equation 

of  state  has  no  effect  on  the  particle  dynamics  defined  by  (3),  or  on  the 

continuum  problem  defined  by  (l)-(2).  For  example,  with  the  polytropic 

equation  of  state  p  =  pQ  +  kp^,  a  measure  of  the  diluteness  of  a  gas 

(in  the  sense  of  kinetic  theory)  is  provided  by  h[u2/k(y  -  l)]1  A'”*1, 

which  estimates  the  average  ratio  of  the  distance  of  closest  approach  of 

* 

two  adjacent  molecules,  to  the  mesh-length. 

For  a  related  similarity  law  involving  mass-spring  models  without  dash- 
pots,  see  Rayleigh,  Proc.  roy.  soc.  A66,  p.  68  (ref.  in  [8,  p.  283]). 
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III.  PLANE  FLOWS 


15.  Elastic  fluids.  The  concept  of  an  elastic  fluid  was  de¬ 
fined  in  I  2.  It  can  also  be  defined  as  an  elastic  continuum  (see  I  18) 
in  which  the  stress  tensor  has  no  shear  component  in  any  direction,  and 
the  pressure  at  any  material  point  is  determined  by  the  specific  volume. 
In  fluid  dynamics,  these  correspond  to  the  assumptions  of  homentropic 
flow.  We  will  treat  below  only  the  case  of  a  homogeneous  elastic  fluid, 
for  simplicity. 

In  terms  of  Lagrangian  independent  variables,  elastic  fluids 
are  defined  by  a  plane  analog  of  Eqs.  (l)-(2).  This  analog  involves  a 
vector  function  x(a;t),  whose  components  are  x(a,b;t)  and  y(a,b;t).  The 
analog  of  (l)  in  this  notation  is 

xa  *b 

This  is  combined  with  a  homentropic  equation  of  state  (2): 

(23)  P  =  PQ  *  F(a)* 

The  function  F(a)  is  chosen  to  agree  with  experimental  values  measured 
under  quasi-static  adiabatic  conditions.  Typically,  one  supposes  that 


(22) 


*tt  ■  -  ‘,Vp’ 


where  a  = 
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x(a;0)  is  given  inside  a  domain  of  interest,  and  the  normal  component 

u  of  x.  =  u  on  the  boundary.  (Alternatively,  one  can  suppose  p  given 
n  ”""t» 

on  a  "free"  boundary.) 

Though  mathematical  existence  and  uniqueness  theorems  seem  not 
to  be  available,  physical  and  numerical  experience  suggests  that  such 
problems  are  usually  well-set,  except  for  the  development  of  shock  waves. 
Shock  waves  must  be  treated  by  special  techniques,  involving  the  Rankine- 
Hugoniot  equations.  These  techniques  will  be  discussed  in  §  17* 

We  will  analyze  below  various  schemes  for  solving  Eqs.  (22) -(23) 
and  related  systems  of  equations  on  digital  computing  machines.  Only  a 
few  of  these  schemes  define  molecular  models  of  matter  in  the  classical 
sense.  But  they  can  all  be  regarded  as  defining  "synthetic  materials," 
or  models  of  matter  which  may  simulate  highly  stressed  real  materials  to 
a  greater  or  less  extent. 

1 6,  Spatial  discretizations.  It  is  easy  to  think  of  spatial 
discretizations  of  the  -system  (22)-(23).  For  example,  one  can  select  a 
square  or  rectangular  network  of  particles  ("molecules")  (a^,  b^)  in 
(a,b)-space,  and  consider  the  particle  paths  x„(t)  =  x(a^,b^;t).  By 
(22) -(23),  these  should  satisfy  the  ordinary  DE’s. 

(24)  x±.  =  o±J  F'(ai(.)  do/dx,  y.j  = 

the  partial  derivatives  being  evaluated  at  (a. ,b  ). 

1  J 

*One  can  also  use  triangular  or  hexagonal  networks  in  place  of  rectangu¬ 
lar  ones. 
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There  are  many  ways  of  evaluating 


xa  % 


ya  yb 


approximately.  For  instance,  one  can  use  the  central  difference  scheme 

(25)  o.  ~  i>j+1  (4AaAb)  . 

1J“  <yi+i,j  -  yi-i,j)(yi,j+i  -h.j-i) 

A  logical  scheme  for  evaluating  da/dx  and  do/dy  is  less  easy  to  select, 
because  of  the  many  options  for  taking  differences  in  (x,y)-space  and 
(a,b)-space.  Unless  care  is  exercised,  one  may  find  that  the  final  for¬ 
mulas  involve  xi±2  ^  and  x±  ^±2,  which  is  undesirable. 

To  avoid  this,  one  can  define  a  nominal  specific  volume  at  each 


half-integer  mesh-point  by 


^Xi,j  “  Xi-1,j^Xi,,3  "  Xi, j-1  ^ 

0i-4’34‘ 


(AaAb)  h 


* 

One  can  then  (for  example)  use  bilinear  interpolations  between  the  four 
members  0  1  .  lt  in  (x,y)-space  or  (a,b)-space,  to  evaluate  Vo  .  Sub- 
stitution  back  into  (24)  gives  a  system  of  second-order  ordinary  DE's  in 
the  vectors  (x. .,y.  ),  in  which  the  right-hand  side  involves  all  eight 

1 J  1 J 

neighboring  mesh-points. 


^Bilinear  interpolation  uses  bilinear  functions  a  +  £x  +  yy  +  &xy  with 
undetermined  coefficients. 


Mechanical  models.  One  can  also  derive  a  spatial  discretization 


by  mechanical  analogy,  placing  a  particle  of  mass  m  at  each  mesh-point, 

joining  adjacent  particles  by  impermeable  but  weightless  straight  rods 

free  to  lengthen  or  shorten  longitudinally,  and  letting  the  "gas"  in 

each  quadrilateral  area  bounded  by  these  rods  expand  and  contract  adia- 

batically.  That  is,  if  N.  ,  (i'  =  i--§>  j'  =  j-^)  is  the  area  of  the 

quadrilateral  with  vertices  at  x(a . ,  .1,  b  i,t),  then  p  ,  ,(t)  = 

1  -2  J  ~2  x  , J 

-  P(N,,  .,/m)  as  in  (23).  One  easily  obtains  formulas  for  the  kinetic 

O  1  y  J 

•  2  *  2 

energy  T  =  Z  m  (x. .  +  y  )/2  and  the  potential  energy 
i j  i j 


V  =  Z  V. ,  . ,  =  -Z  / p.  ,  dA. ,  of  such  a  system.  If  p  =  kp^  = 

1  /  J  1  j  J  1  j  J 

P0(Aq/A)^  (7  ^  l),  then  /pdA  =  kA1  with  k  =  pqA£,  ior  example. 

As  in  Lagrangian  dynamical  systems  without  external  forces  generally, 
the  equations  of  motion  are  then  given  variationally,  by  setting 
d(dL/d4j)/dt  =  dl/dq^,  where  L  =  T  -  V.  This  model  was  suggested  by  one 
of  the  authors  (Garrett  Birkhoff,  unpublished  note,  1959). 

A  program  for  computing  the  motion  of  the  preceding  "synthetic 
fluid,"  with  finite  time- step s> has  been  devised  by  Walter  Goad  [25].  His 
method  differs  in  several  respects  from  an  approximate  method  for  a 
slightly  different  model  suggested  earlier  by  Harwood  Kolsky  [14].  It  is 
not  clear  that,  in  the  limit  Atl  0,  Kolsky’ s  model  corresponds  to  a 
Lagrangian  system.  Lagrangian  systems  are  not  only  esthetically  satisfy¬ 
ing,  they  have  the  advantage  that  they  make  L(q,q)  constant  in  time, 


*  _  . 

For  the  continuum  problem  of  s  14,  similar  variational  formulations 
have  been  given  by  H.  Bateman,  "Partial  differential  equations," 
p.  l64;  Lamb,  p.  188;  P.  Lush  and  T.  Cherry  QJMAM  9  (1956),  6-21. 


thus  providing  a  useful  check  on  computations. 


17.  Shocks  and  viscosity.  The  primary  aim  of  the  computation 
schemes  devised  by  Kolsky  and  Goad  was  to  calculate  approximately  the 
motion  of  an  elastic  fluid,  allowing  for  shock  discontinuities.  Much  as 
in  I  6,  one  does  this  by  introducing  a  large  so-called  artificial  vis¬ 
cosity,  chosen  so  as  to  smooth  out  irregularities  behind  the  shock  while 
permitting  a  step  shock  front  to  propagate  without  change  of  form.  We 
yiH  now  try  to  correlate  such  schemes  with  the  concept  of  a  viscous 
elastic  fluid. 

A  viscous  elastic  fluid  may  be  defined  as  a  continuum  in  which  the 
stress  tensor  (matrix)  | | |  is  a  single-valued  linear  function  of  the 
rate-of- strain  tensor  | l^u^/dx^l |  at  any  point.  We  consider  below  only 
homogeneous  and  isotropic  viscous  elastic  fluids,  which  necessarily 
satisfy  the  Navier-Stokes  equations. 

Viscous  effects  in  such  fluids  are  known  to  be  characterized  by 
their  shear  viscosity  p  and  their  bulk  viscosity  p*,  material  constants 
which  depend  on  the  temperature  T  and  pressure  p  (the  "state"  (p,T))  of 
the  fluid.  To  determine  exactly  the  motion  of  a  viscous  elastic  fluid, 
one  must  know  its  conductivity  k  =  k(p,T)  as  well  as  p  and  P*.  In  prac¬ 
tice,  one  is  satisfied  with  approximate  results,  and  makes  simplifying 

assumptions  such  as  p  =  const.,  p'  =0. 

In  compressible  flow  computations,  one  seldom  attempts  to  simulate 
real  physical  viscosity  at  all.  Instead,  since  the  assumption  p  >  0 
would  give  rise  to  the  complication  of  no  slip  on  solid  boundaries,  one 


tends  rather  to  set  p  =  0  and  to  introduce  a  large  artificial  bulk 


viscosity  p*»  p  in  calculating  plane  flows.  Bulk  viscosity  also  has 
the  advantage  over  shear  viscosity  of  introducing  a  scalar  pressure 
stress,  and  not  a  symmetric  tensor  shear  stress  with  three  components, 
in  plane  flows.  The  whole  purpose  of  p*  is  to  smooth  out  "random"  fluc¬ 
tuations  behind  shocks,  and  to  make  the  calculated  "shock  thickness"  a 
few  mesh- lengths. 

As  stated,  the  real  aim  is  thus  to  simulate  an  elastic  fluid 
with  shocks,  in  which  the  shocks  satisfy  Rankine-Hugoniot  relations  de¬ 
termined  by  a  function  p(l,a),  expressing  the  pressure  at  any  point  as 
a  function  of  the  internal  energy  per  unit  mass  of  the  fluid,  I,  and  the 
local  specific  volume.  In  a  perfect  gas,  I  =  T  and  p  =  RT/cr  =  Rl/c^cr. 

In  practice,  this  simulation  has  been  successfully  achieved  in 
many  cases,  using  finite  time-steps  At  governed  by  variants  of  the 
Courant  stability  condition  [14,  p.  1 3] .  It  has  the  advantage  over 
"shock  fitting"  that  it  treats  all  points  alike.  It  has  the  disadvan¬ 
tage  of  requiring  a  second  dependent  variable  (internal  energy  I  or  en¬ 
tropy  S)  to  define  the  local  "state"  of  the  fluid. 

18.  The  Kolsky-Goad  fluid.  To  apply  the  computational  methods  of 
Kolsky  and  Goad  to  real  materials,  one  must  know  both  the  pressure 
p(a,T)  and  the  internal  energy  l(cr,T)  as  functions  of  the  variables  cr,T. 
Especially,  one  must  know  p(cr,l)  across  a  single  shock  (the  Hugoniot 

Real  bulk  viscosity  is  subject  to  many  anomalies;  only  in  monatomic  gas 
is  Stokes'  assumption  p'  =  0  justified. 


curve  of  [26,  p.  152]).  For  solids,  this  will  depend  on  the  effective 
Poisson  ratio,  which  may  be  unknown  under  the  conditions  considered;  see 
§  19.  This  uncertainty  suggests  a  number  of  interesting  questions. 

It  is  also  interesting  to  consider  the  "Kolsky-Goad  fluid"  defined 
by  the  mechanical  model  of  §  16  as  a  synthetic  material.  For  it  to  be 
stable,  p(a)  must  be  a  decreasing  function. 

Though  the  microscopic  appearance  of  this  spatial  discretization 
of  ( 22 ) — ( 23 )  is  quite  unlike  that  of  real  liquids  and  gases,  its  macro¬ 
scopic  behavior  is  that  of  an  elastic  fluid.  One  wonders  whether  it 
might  not  be  used,  for  example,  to  simulate  real  equations  of  state  for 

gases  and  liquids,  including  phase-transition.  For  instance,  if  one  de- 

~~2 

fines  the  temperature  T  as  the  mean  kinetic  energy  mv  /2  of  translation 
of  the  particles,  one  can  ask  for  p(o,T). 

Due  to  the  complications  which  arise  when  the  quadrilateral  areas 
A.  become  extremely  distorted  (which  are  also  inconvenient  in  ordinary 

J.  J 

calculations),  this  may  not  be  easy.  But  one  wonders  whether  some  modi¬ 
fied  scheme,  perhaps  using  triangular  or  hexagonal  areas  and  allowing 
unwanted  rods  to  disappear  and  new  rods  to  subdivide  irregular  areas, 

might  not  be  successful  in  simulating  the  behavior  of  real  gases  and 

* 

liquids.  It  must  be  confessed  that  other  models  for  liquids  also  fail 
to  simulate  real  molecular  physics  in  any  detail. 

*E.g.,  the  "free  volume"  and  "hole"  theories  [17,  §  4.4,  1  4.8].  See 
also  J.  Frenkel,  "Kinetic  theory  of  liquids, "Oxford,  1946;  M.  Born  and 
H.  S.  Green,  "A  general  kinetic  theory  of  liquids,"  Cambridge  Univ. 
Press,  1946. 


19.  Elastic  solids.  In  continuum  mechanics,  one  defines  an  elas¬ 


tic  solid  as  a  continuum  in  which  the  stress  tensor  | |p^j | |  is  a  single¬ 
valued  function  of  the  strain  tensor  | | bx^/ba  | | .  The  theory  of  elastic 
solids  is  highly  developed  only  for  infinitesimal  deformations  from  sta¬ 
tic  equilibrium  of  homogeneous ,  isotropic  media  which  have  not  been  pre¬ 
stressed.  One  then  assumes  |  |p^j|  |  =  0  if  bx^/ba^  =  5^,  and  makes  |  |p^|  I 
depend  linearly  on  the  difference  | Idx^/da^ | |  -  | |6^j | | . 

Under  these  assumptions,  the  elastic  characteristics  of  a  three- 
dimensional  medium  can  be  shown  [21 ]  to  be  characterized  by  a  rigidity 
or  shear  modulus  p  and  a  bulk  modulus  p’  =  A  +  (2p/3 );  the  analogy  with 
the  viscosity  coefficients  is  obvious.  In  the  plane,  one  has  p1  =  A  +  p 
instead. 


It  is  frequently  suggested  that  the  equations  of  sin  elastic  fluid 
can  be  substituted  for  those  of  elastic  solids,  in  two  cases:  (i)  in 
one-dimensional  motion  [5,  i§  3-5  and  li  97-101],  and  (ii)  in  violent  ex¬ 
plosions,  when  stresses  far  exceed  the  yield  point  of  the  solid  considered. 
However,  this  substitution  must  be  employed  with  caution,  because  of  am¬ 
biguity  as  to  the  relevant  modulus  of  elasticity. 

Thus,  across  a  slab,  the  velocity  of  compressive  acoustic  waves 

JL  JL. 

is  [(A  +  2p)/p]2;  along  a  bar,  it  is  [E/p]2  according  to  the  simple 
theory,  where  [21,  p.  42 ]  E  =  p  [(3A  +  2p ) / ( A  +  p)].  For  accurate  results, 
edge  corrections  are  needed  [26,  p.  35].  Likewise  the  triaxial  stress 
required  to  produce  a  uniaxial  strain  (e,0,0)  is  (A  +  2 p,  A,A)e,  and  this 
is  very  different  from  the  hydrostatic  pressure  (A  +  2p/3)e  required  to 
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produce  the  same  "bulk,  compression.  Thus  for  steel,  in  which  p/4  =  1l/l4 
because  the  Poisson  ratio  v  =  \/ 2(A  +  p)  =  0.28  [26,  p.15^1*  x“ 
component  of  pressure,  ^-s  nearly  twice  as  great  in  the  first  case 

as  in  the  second. 

All  this  shows  that  elastic  solids  are  much  more  complicated  than 
fluids  in  their  reactions  to  one-dimensional  stresses. 

As  regards  waves  of  finite  amplitude,  consider  a  shock  front. 
Lateral  contraction  over  the  short  transition-time  would  produce  unbe¬ 
lievably  large  acceleration  stresses;  hence  the  compress ion- jump  is  al- 

* 

most  surely  uniaxial.  In  view  of  Bridgman's  data  on  the  increase  in 
shear  strength  under  hydrostatic  pressure  and  under  rate  of  strain,  how 
can  one  be  sure  that  shear  stresses  (proportional  to  p  in  the  elastic 
range)  are  negligible? 

Because  of  all  these  uncertainties,  it  seems  desirable  to  ex¬ 
plore  the  many  kinds  of  matter  which  have  been  found  in  Nature,  and 
molecular  models  for  these  kinds  of  matter.  Besides  elastic  fluids  and 
solids,  we  have  already  mentioned  viscoelastic  fluids  and  Maxwell  solids. 
Viscoelastic  and  plastic  solids  should  certainly  also  be  considered. 
Perhaps  one  should  also  consider  glasslike,  rubberlike,  and  doughlxke 
substances.**  What  about  sand  and  sintered  metals? 

*P.  J.  Bridgman,  "Studies  in  large  plastic  flow  and  fracture,"  McGraw- 
Hill,  1952,  esp.  Ch.  16  (see  also  Ch.  12,  where  various  complicated 
phenomena  affecting  volume  changes  are  described) . 

**We  follow  the  classification  of  R.  Houwink,  "Elasticity,  plasticity, 
and  structure  of  matter,"  2d  ed.,  Harren  press,  1953.  Houwink  also 
points  out  distinctive  mechanical  properties  of  casein,  starch,  and 
clay.  Sols,  gels,  paints,  etc.,  have  other  peculiarities. 
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Against  the  background  of  this  general  philosophy,  we  will  now  ex¬ 
plore  a  few  very  simple  mass-spring  models,  which  seem  especially  promis¬ 
ing  as  possible  substitutes  for  pure  fluids  in  computation  schemes  which 
purport  to  approximate  the  motion  of  violently  stressed  solids. 

20.  Mass-spring  networks.  Mass-spring  networks  in  the  plane  and 
in  space  have  been  studied  since  the  time  of  Cauchy,  as  providing  mole¬ 
cular  models  of  elastic  solids.  In  accordance  with  the  terminology  in¬ 
troduced  earlier,  one  can  think  of  them  as  defining  synthetic  crystals. 

It  is  interesting  to  consider  the  mechanical  properties  of  such 
synthetic  crystals  for  their  own  sake,  especially  when  one  really  wishes 
to  simulate  a  solid.  Following  Cauchy  and  Kelvin,  Brillouin  [3]  has 
investigated  the  propagation  of  waves  through  such  systems  connected  by 
1 i nsar  springs  (7  =  -1 ) ;  however,  for  hydrodynamical  computations,  it 
seems  more  appropriate  to  choose  F(cr)  to  match  p  in  (23)  for  uniform  (hy¬ 
drostatic)  compression.  For  example,  for  F(cr)  =  ka  ^ ,  one  might  consider 
the  non-linear  force-law  f(r)  =  -a  r"  .  A  simple  dimensional  analysis 
shows  that  this  makes  B  =  2y  -  1 ;  the  relation  between  k  and  a  depends 
on  the  type  of  crystal. 

Simple  networks.  The  simplest  plane  crystals  have  polygonal  cells; 
the  cases  of  squares,  equilateral  triangles,  and  regular  hexagons  are ( 
the  most  interesting.  If ,  in  a  square  (or  hexagonal) 
network,  each  particle  acts  only  on  its  4  (or  3) 
immediate  neighbors  as  indicated  in  Fig.  1 ,  then 
the  network  is  obviously  unstable  in  compression. 


because  the  area  enclosed  can  be  flattened  without  shortening  the 
springs.  Moreover,  the  equation  of  state  for  this  model  is  not  a  func¬ 
tion  of  specific  volume  —  if  the  area  of  a  cell  goes  to  zero  without 
shortening  the  springs,  the  pressure  does  not  change. 

Though  unstable  in  compression,  the  square  network  defined  above 
is  stable  in  tension.  This  shows  that  the  dynamic  behavior  of  plane 
mass-spring  systems  may  change  when  pQ  is  changed  in  the  equation  of 
state  simulated  —  even  though  this  is  not  true  of  continuous  fluids  or 
for  networks  in  one  space  dimension. 

Braced  square  networks.  Another  synthetic  crystal  is  defined  by 
the  diagonally  braced  square  cell  of  Fig.  2.  Each  net  point  is  connected 
with  its  8  nearest  neighbors.  Though  not  necessarily  isotropic,  such 
crystals  have  two  orthogonal  "principal  planes"  of  reflection  symmetry. 

To  compute  the  elastic  constants  of  this 
model,  we  suppose  that  the  spring  con¬ 
stants  for  the  horizontal  and  vertical 
springs  are  a,  while  those  for  the  diag¬ 
onal  springs  are  b.  The  associated 
stiffness  constants  are  then  a  and  b//2,  respectively. 

The  elastic  moduli  can  be  easily  computed.  We  obtain  c1 1  =  a  +  b, 

c  =  b.  and  =  b.  The  stress-strain  relation  is  given  by 
12  66 


*We  use  the  notation  of  I.  Sokolnikoff,  "Theory  of  elasticity,"  second 
edition,  McGraw-Hill,  1956. 
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where  t. .  and  e#J  are  components  of  stress  and  strain.  The  condition  for 
isotropy  is  c1 1  =  c  +  2Cgg;  a  mass-spring  system  of  this  type,  in  which 
a  =  2b,  is  elastically  isotropic.  The  Poisson  ratio  v  =  b/(a  +  b)  is  l/? 
in  the  isotropic  case  as  it  should  be,  for  reasons  which  will  now  be  ex¬ 
plained. 


21 .  Networks  as  spatial  discretizations.  In  one  dimension,  mass¬ 
spring  systems  define  useful  spatial  discretizations  of  elastic  solids 
(and  fluids).  It  is  natural  to  ask  if  this  is  also  true  in  the  plane 
(and  in  space). 

As  far  as  infinitesimal  deformations  of  homogeneous  isotropic  solids 
are  concerned,  this  question  was  partly  answered  by  Cauchy.  He  assumed 
implicitly  what  may  be  called  Cauchy’s  Hypothesis:  that  the  mass-spring 
system  remained  in  static  equilibrium  under  any  affine  deformation.  This 
is  true,  with  central  forces,  for  any  crystal  having  central  symmetry  in 
all  its  mass-points  —  hence  for  all  the  models  of  I  20. 

Under  Cauchy's  hypothesis,  p.  =  A,  which  implies  v  =  1/3  in  the  plane 
(as  above),  and  v  =  1/4  in  space.  But  experimentally,  v  is  fair  from  l/4 
in  many  solids;  thus  for  lead,  v  =  0.45  [2 6,  p.  154]. 

It  follows  that  one  cannot  approximate  the  behavior  of  most  elastic 
solids  by  simple  mass-spring  networks  centrally  symmetric  in  every  mass- 
point,  even  as  far  as  plane  waves  of  infinitesimal  amplitude  are  con¬ 
cerned.  In  particular,  the  bulk  modulus  and  shear  (rigidity)  modulus 
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are  proportional,  so  that  one  cannot  easily  simulate  an  elastic  fluid. 
This  is  also  true  of  the  ionic  crystals  (with  oppositely  charged  par¬ 
ticles)  used  to  simulate  the  alkali  halides  ([1],  [2],  [3])  —  though 
various  interesting  properties  of  real  crystals  ( e.g. ,  "Brillouin  zones") 
are  simulated  qualitatively. 

Anisotropy.  The  mechanical  behavior  of  synthetic  crystals  under 
small  deformations  is  determined  by  their  moduli  of  elasticity.  Hence, 
the  braced  square  network  discussed  earlier  is  elastically  isotropic, 
even  though  no  molecular  model  can  be  isotropic  in  all  its  physical  pro¬ 
perties.  However,  even  elastic  isotropy  is  lost  under  large  defor¬ 
mations,  a  fact  which  makes  it  difficult  to  use  synthetic  crystals  to 
simulate  isotropic  polycrystalline  solids  (like  steel)  with  v  near  1  /4. 

Consider  for  example  a  simple  synthetic  crystal  corresponding  to 

a  mass-spring  network  of  equilateral  triangles.  Let  the  springs  satisfy 

-6  “7 

the  force-law  f(r)  =  a  r  .  The  pressure  will  then  be  p  =  k  a  if 

6  =  2y  -  1  and  if  a  =  k(^/4m)-^//^  where  m  is'  the  mass  of  a  "molecule." 

Then, under  large  uniaxial  compressions,  one  obtains  the  configura¬ 
tion  of  Fig.  3.  It  is  easily  verified 
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Also  this  crystal  is  isotropic  for  infinitesimal  deformations;  it 
is  highly  non-isotropic  (and  unsuitable  for  hydrodynamical  computations) 
as  regards  finite  deformations. 

Kinetic  theory.  As  in  the  one-dimensional  case,  it  would  be  in¬ 
teresting  to  compute  equations  of  state  for  synthetic  crystals  subject  to 
various  force-laws,  expressing  the  specific  volume  0  =  a(p,T)  as  a  func¬ 
tion  of  the  temperature .  As  in  §  10,  one  would  have  to  define  T  in  terms 
of  statistical  equilibrium  with  a  perfect  gas  at  the  stated  temperature, 
and  not  in  terms  of  the  mean  energy  per  particle.  Expansion  anisotropy 
could  also  be  computed. 

22.  Universal  force-laws.  In  the  computation  schemes  and  syn¬ 
thetic  crystals  discussed  so  far,  forces  have  been  assumed  to  act  only 
between  preassigned  "near  neighbors".  It  is,  of  course,  more  reasonable 
to  let  forces  be  determined  entirely  by  the  distances  between  the 
particle-pairs  involved,  irrespective  of  their  initial  locations.  (For 
the  small  deformations  of  classical  elasticity  theory,  to  be  sure,  the 
distinction  is  largely  academic.) 

The  artificial  restriction  of  molecular  forces  to  preassigned 
pairs  of  neighbors  has  been  avoided  by  Born  ([1],  [2])  using  Made lung's 
constants.  However,  he  makes  the  equally  artificial  assumption  of  per¬ 
fectly  periodic  displacements  in  calculating  the  mechanical  properties 
of  crystals.  (Somewhat  similarly,  Cauchy  and  others  have  made  calcula¬ 
tions  assuming  that  all  particles  except  the  one  under  consideration  re¬ 
mained  stationary.) 
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Only  recently  have  genuine  "kinetic  theory"  calculations  been  made, 

in  which  forces  were  allowed  to  act  between  arbitrary  pairs  of  suffi- 

* 

ciently  near  particles,  by  F.  Harlow  and  (for  the  model  discussed  in 
I  1 6)  by  Harwood  Kolsky. 

However,  Harlow's  model  uses  central  forces  and  a  simple  initial 
configuration;  it  seems  likely  that  it  has  shear  rigidity  if  it  is  iso¬ 
tropic.  Moreover  it  retains  the  artificiality  of  not  having  a  universal 
force-law,  applicable  to  all  sufficiently  close  particle -pairs.  It 
would  seem  to  us  that  the  device  of  having  a  cutoff  radius  (depending 
on  the  entropies  S  or  temperatures  T  of  the  particle-pairs  involved) 
would  be  more  logical  and  perhaps  as  effective  computationally  as  re¬ 
stricting  forces  to  a  preassigned  number  of  neighbors  of  each  particle. 
This  could  be  combined  with  an  artificial  conductivity,  acting  between 
"nearby"  particle-pairs  with  f >  0,  so  as  to  smooth  out  both  shocks 

1 J 

and  irregularities  due  to  particle  diffusion. 

Postulating  forces  with  such  a  cutoff  radius  was  originally  sug¬ 
gested  by  Newton  (loc.  cit.  in  §  3).  It  avoids  the  possible  paradox  of 
infinite  potential  energy  associated  with  distant  particles,  implicit  in 
some  universal  force-laws.  For  instance,  if  f(r)  =  kr  p  with  p  <  4, 
then  the  volume  integral  C/rP-1  r2dr(C  =  4jck/(p  -  1))  of  probable  poten¬ 
tial  energy  of  randomly  distributed  particles  is  infinite.  This  paradox 

*F.  H.  Harlow  and  B.  D.  Meixner,  "The  particle-and-force  computing 
method  for  fluid  dynamics,"  Report  LAMS-2567>  Bos  Alamos,  19^1 . 

**H.  Kolsky,  "The  nearest  neighbor  hydrodynamics  calculation,"  Memo, 
dated  July  1 ,  1 9^1 . 


-52- 


has  been  discussed  by  Maxwell  (see  [8,  1  1 65] )• 

In  the  computation  schemes  of  Harlow  and  Kolsky  mentioned  above, 
an  artificial  viscosity  is  introduced,  and  a  time-dependent  entropy 
is  attributed  to  particles  (or  cells).  It  would  be  interesting  to  study 
the  synthetic  materials  defined  by  the  same  models  without  viscosity, 
and  to  compute  the  equations  of  state  p(o,T)  and  internal  energy  functions 
l(c,T)  for  the  synthetic  materials  which  they  define.  This  can  be  re¬ 
garded  as  an  extension  to  the  plane  and  to  universal  force-laws  of  the 
ideas  introduced  in  i  11,  for  one- dimensional  mass-spring  systems. 

Plastic  solids.  Under  increasing  shear  strains,  particle-and- 
force  systems  with  universal  force-laws  behave  somewhat  like  plastic 
solids.  Consider  for  example  the  square  network  of  1  20  with  molecular 
spacing  h  and,  for  simplicity,  a  universal  force-law  with  cutoff  dis¬ 
tance  d  =  2h.  For  small  deformations  from  equilibrium,  the  system  will 
behave  like  a  braced  square  network.  But  if  a  shear  strain  is  imposed 
which  slides  each  m-th  row  through  a  distance  mh/2,  then  symmetry  shows 
that  there  will  be  static  equilibrium  (zero  stress). 

It  would  be  interesting  to  try  to  simulate  other  features  of  plas¬ 
tic  behavior  by  quasi-molecular  particle-and-force  systems.  Especially, 

it  would  be  interesting  to  determine  the  relations  to  large  strains  of 

* 

ideal  simple  crystals  corresponding  to  the  fourteen  Bravais  lattices. 

•X- 

C.  K.  Hel,  "introduction  to  solid  state  physics,"  Wiley,  1956,  Ch.  I. 

To  simulate  the  effect  of  dislocations  in  imperfect  crystals  would 
clearly  be  harder.  See  R.  Hill,  "The  mathematical  theory  of  plasticity," 
Oxford,  1956,  Ch.  I. 
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23.  PIC  models.  Perhaps  the  most  novel  scheme  which  has  been 


used  successfully  to  simulate  an  elastic  fluid  with  shocks  is  the 

particle-in-cell  or  PIC  computation  scheme  developed  at  Los  Alamos  by 

Harlow  and  his  co-workers.  As  the  usefulness  of  this  scheme  as  a  compu- 

* 

tational  tool  has  been  discussed  in  some  detail  in  other  reports,  we 

will  consider  here  primarily  the  synthetic  materials  which  its  variants 

define.  (The  scheme  now  used  in  hydrodynamical  computations  evolved 

** 

after  extensive  trial-and-error  experiments  with  such  variants.  ) 

To  represent  a  homogeneous  material,  a  box  is  supposed  divided 
into  square  (or  rectangular)  cells,  each  of  which  contains  a  number 
Nc(t)  of  identical  particles.  Each  cell  C  has  density  pc  =  Nc  (in  suit¬ 
able  units)  and  pressure  p  =  f(N  ,1  ),  where  I  is  an  internal  energy 

c  c  c  c 

depending  on  the  past  history  of  particles  which  have  been  in  the  cell. 

In  the  simple  adiabatic  case,  we  would  have  p  =  f(N  ),  but  the  various 

c  c 

averaging  processes  used  in  actual  PIC  codes  (KAREN,  SUNBEAM,  etc.)  act 
like  artificial  viscosity  in  converting  mechanical  energy  into  heat. 

An  obvious  peculiarity  of  the  synthetic  material  defined  by  some 
such  schemes  is  the  instability  of  static  equilibrium  unless  N  is  a 
half-integer.  Unless  the  cells  neighboring  a  given  cell  on  each  side 
contain  the  same  number  of  particles,  the  cell  will  be  accelerated.  It 
follows  that  in  any  stable  static  array  of  particles,  all  "red"  squares 

*See  [13],  [24],  and  F.  H.  Harlow,  J.  assoc,  comp,  mach.4  ( 1 957) >  137-42. 
Described  in  earlier  Los  Alamos  Reports  LAMS- 1 95°  and  LAMS-2082. 


-54- 


of  a  checkerboard  must  contain  the  same  number  of  particles,  and  so  must 
all  "black"  squares.  The  converse  is  also  true,  which  shows  that  static 
stability  is  possible  if  and  only  if  the  average  density  p  =  Nq  is  an 
integer  or  half-integer. 

Obviously,  the  "diophantine "  static  instability  described  above  be¬ 
comes  less  serious  as  increases.  It  is  interesting  to  consider  the 

behavior  of  PIC  material  in  the  limit  N  -*  ».  In  one  dimension  this  can 

c 

be  treated  semi-analytically,  taking  the  density  p(x,t)  as  dependent  var¬ 
iable.  Possibly  convergence  to  the  behavior  of  an  elastic  fluid  can  be 
proved  in  this  way,  taking  the  iterated  limit  as  Nc  and  cell  size 
tends  to  zero. 

It  has  been  shown  [13,  p.  16]  that,  in  one  dimension,  the  PIC  ma¬ 
terial  has  an  effective  (bulk)  viscosity  p*  =  p|p|s/2  (s  =  cell  length); 
the  analogy  with  the  formula  p  =  pcu/3  of  the  kinetic  theory  is  obvious. 
On  the  other  hand,  the  natural  shear  viscosity  is  zero:  the  drag  layer 
along  slip  boundaries  is  at  most  one  cell  wide,  even  if  averaging  is 
used. 

The  PIC  material  is  highly  non-linear  (because  of  its  "diophantine" 
behavior,  remotely  suggestive  of  quantum  mechanics).  Hence  it  would  be 
interesting  to  compare  the  formula  for  bulk  viscosity  quoted  above  with 
attenuation  rates  for  sinusoidal  waves  of  variable  frequency  and  ampli¬ 
tude,  in  directions  parallel  to  cell  boundaries  and  at  a  k5°  angle  with 
them. 

*F.  Harlow,  J.  assoc,  comp,  mach.4  (l957)>  P*  139* 
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24.  Synthetic  materials  testing.  As  the  number  of  "synthetic 
materials"  proposed  for  hydrodynamic  computations  increases,  the  need  for 
objective  standards  of  comparison  clearly  grows.  It  seems  to  us  that 
standard  materials  testing  procedures  used  in  mechanical  engineering  pro¬ 
vide  good  precedents  in  this  connection.  (Only  dynamic  tests  would  ordi¬ 
narily  be  of  interest.) 

Thus,  in  the  linear  range,  one  can  derive  moduli  of  elasticity  and 
viscosity  by  studying  the  velocity  of  propagation  and  the  rate  of  atten¬ 
uation  of  plane  waves  in  various  directions.  Explosions  of  highly  com¬ 
pressed  spheres,  cubes  (squares),  etc.,  give  some  useful  comparisons  in 
the  non-linear  range.  So  do  simple  plane  shocks,  of  varying  intensity 
and  moving  with  varying  velocity.  Impact  tests  (say  of  high-speed 
rigid  missiles),  somewhat  analogous  to  drop  tests,  could  also  be  used. 

But  just  as  in  the  case  of  ordinary  materials  testing,  the  stan¬ 
dardization  of  such  a  series  of  tests  can  probably  best  be  decided  by  a 
committee.  Therefore  we  will  not  elaborate  further  on  the  subject. 


*Though  the  problems  which  they  attempt  to  solve  are  invariant  under 
transformation  to  moving  axes,  the  same  is  not  true  of  digital  Eulerian 
or  quasi-Eulerian  schemes,  such  as  PIC. 
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